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A.  Statement  of  Problem  Studied 

Contract  No.  MIPR-ARO  124-84  was  a  $35,000  grant  to  investigate  methods 
for  improving  one  and  two  dimensional  shock  flow  calculations  by  combining 
several  approaches  such  as  mesh  refinement  and  moving  grids  for  optimal 
accuracy  and  efficiency. 

B.  Summary  of  Important  Results 

Extremely  accurate  shock  flow  results  were  obtained  for  standard  one 
dimensional  model  problems.  Inviscid  hydrodynamics  are  piecewise  continuous; 
derivatives  at  a  discontinuity  must  be  formed  by  appropriate  one-sided 
difference  approximations  consistent  with_calculus.  Because  the  spatial 
derivative  at  a  discontinuity  is  unbounded,  it  can  be  made  finite  if  the 
discontinuity  is  spread  out  using  a  consistent  step  function  approximation, 
for  example,  an  arctan.  The  derivative  with  the  discontinuity  is  a  finite 
approximation  to  a  delta-function. 

Second,  a  moving  grid  scheme  was  used.  The  grid  velocity  which  is 
arbitrary  is  taken  to  be  that  velocity  in  which  the  set  of  dependent  variables 
appear  steady  in  the  least  squares  sense.  At  a  shock,  contact  surface  or 
rarefaction,  the  grid  velocity  for  each  component  was  equal,  but  not 
necessarily  equal  in  each  continuous  subregion.  Wave  steepening  drives  mesh 
points  together  and  whenever  the  distance  of  closest  approach  fell  below 
0.0001  and  the  gradients  in  that  region  grew  from  successive  time  steps,  a  new 
discontinuity  was  inserted.  A  conservative  static  rezoning  scheme  was  used  as 
a  post  processor  after  every  time  step  permitting  grids  to  move  without 
constraints. 

After  extensive  testing,  it  was  found  that  the  standard  local  three  finite 
difference  scheme  performed  very  poorly  for  derivative  estimates  in  the 
rarefaction  region.  However,  when  a  global  monotonic  spline  interpolation 
developed  by  Fritsch  and  Carlson  was  used,  vastly  improved  accuracy  was 
obtained. 


Mesh  refinement  was  used  for  interacting  shock  waves  in  conjunction  with  a 
Riemann  solver.  The  Riemann  solver  was  used  only  to  predict  the  outcome  of  a 
wave  interaction,  not  to  advance  the  solution  in  time.  Mesh  points  were 
dynamically  added  to  the  computational  domain  for  new  shocks,  contacts,  and 
rarefaction  waves.  It  was  thus  possible  to  follow  the  interaction  of  shocks 
with  both  step  up  and  step-down  contact  interactions  using  a  minimum  number  of 
points  extremely  accurately. 

Finally  a  two  dimensional  freely  propagating  shock  was  calculated  using 
the  basic  concepts  developed  in  one-dimension.  In  20,  the  components  of  the 
grid  velocity  were  found  by  rotation  into  the  direction  of  the  maximum 
gradient,  and  a  translation  to  a  frame  in  which  the  dependent  variables  are 
steady.  Appropriate  one  sided  differences  were  used  at  the  shock.  When 
standard  stencils  could  not  be  used  at  or  near  the  shock,  a  second  order 
Taylor  series  on  scattered  data  was  used.  The  results  for  this  problem  were 
excellent. 

Because  of  the  limited  contract  grant,  no  mathematical  analysis  was  done. 
The  extension  to  the  general  case  of  20  interacting  shocks  is  unfinished,  but 
the  essential  building  blocks  have  been  obtained.  A  superior  method  for 
differencing  in  two  or  three  dimensions  irregardless  of  a  grid  data  structure 
appears  to  be  Hardy's  [1]  multiquadric  interpolation  (MQI)  scheme. 

MQI  is  a  locally  global  interpolation  scheme  which  uses  upper  hyperboloids 
as  basic  functions.  Franke  [2]  states  that  of  all  the  methods  he  has  tested 
MQI  is  outstanding  in  accuracy,  visual  appearance,  ease  of  use,  and  timing. 

MQI  is  continuously  differentiable  permitting  arbitrary  data  location  without 
the  dispersion  and  dissipation  truncation  errors  associated  with  local  low 
order  interpolation  based  difference  schemes.  MQI  has  the  power  and  accuracy 
of  spectral  methods  without  the  restriction  of  a  logically  rectangular  tensor 
product  formulation. 

Whenever  possible,  the  analysis  of  Ben-Dor  and  Glass  [3]  and  Mirels  [4] 
will  be  used  for  the  2D  regular  and  Mach  reflection  problem.  However,  the 
synthesis  of  the  various  aspects  has  yet  to  be  done,  and  will  be  postponed 
until  future  funding. 
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Aba tract 


Five  one-diaan«ional  (ID)  and  one  2D  ahock  wave 
problems  which  propagate  obliquely  to  the  coordinate  axea 
are  aolvad  by  a  second-order  tine -marching  method.  The 
aolution  region  ia  aaauaed  to  be  piaceviae  continuoua, 
with  any  "discontinuities"  which  may  develop  being 
rapraaanted  by  an  arccan  approximation  to  a  atap  func¬ 
tion.  fnmaiiiataly  behind  or  ahead  of  a  flagged  "discon¬ 
tinuity",  appropriate  one -aided  derivativea  are  uaad. 

An  explicit  moving  grid  technique  ia  combined  with  the 
tine-incegracion  achene  which  yielda  the  correct  velocity 
at  "discontinuities".  Shock  fitting  ia  noc  only  aimply 
handled,  but  it  ia  handled  automatically  and  correctly 
by  the  choice  of  the  grid  velocitiea.  The  regularixation 
problem  aaaociated  with  moving  grida  ia  handled  by  a  re- 
toning  baaed  on  equidiatributing  the  component  averaged 
third  derivative.  Becauae  the  approximate  atep-function 
and  the  above-mentioned  ruled  for  derivative  formation 
are  mathematically  conaiatent,  artificial  viaeoaity  ia 
unneceaaary.  The  dame  aecond-order  time  integration 
acheme  ia  uaed  throughout  the  entire  apatial  domain.  For 
accurace  phyaica  of  wave  interactional  a  nonlinear 
Riemann  solver  ia  applied.  In  conjunction  with  Che 
Riemann  aolvcr,  adaptive  neah  refinement  injecta  maah 
points  into  the  computational  domain  automatically  to 
define  a  poat  colliaion  structure. 

The  Tine-Marching  Scheme 

For  brevity,  the  reader  ia  referred  to  Richtmayer  and 
Morton  [1]  for  deeaila.  The  two-dimensional  conservation 
law  form  of  the  fluid-dynamic  equationa  ia 

He  *  lx  ♦  £y  “  °,  (1) 

where  the  subscripts  t,  x,  and  y  denote  partial  deriva¬ 
tives.  The  fluxes  F  and  C  are  related  to  the  conserva¬ 
tive  variables  U  in  the  following  manner: 

lx  *  *  Hx!  £y  ■  B  Uy  (2) 

where  A  and  8  arc  the  appropriate  Jacobians. 

However,  the  matrix  representation  of  the  Jacobiana 
can  be  ambiguous  regarding  Che  nonlinear  matrix  elements. 
The  following  procedure  was  uaed  to  guarantee  rigorous 
equalicy.  Given  F*.  U*.  Gy  and  Uy  aa  known  quantities, 
the  unknown  nonlinear  matrix  elements  of  the  Jacobiana, 

A  and  3  are  solved. 

The  solution  U,  at  the  new  time  t°*^  is  given  by  a 
second  order  Taylor  aeries  expansion  aa 


U(n-l)  -  U(n)  -  /  dt  {£♦  C  } 

-x  -y 

♦  \  //dt  dt*  {[AfF^  ♦  ♦  (Bl^  ♦  C^Iy) 


(3) 


*  Work  performed  under  Che  auspicea  of  the  U.S.  Depart¬ 
ment  of  Energy  by  the  Lawrence  Livermore  National  Labora¬ 
tory  under  contract  No.  W-7405-ENC-48  and  partially 
supported  by  the  Army  Research  Office  contract  No. 
MIPR-AKO  124-84. 


The  Moving  grid  Scheme 

The  approach  taken  in  this  paper  is  a  modification 
of  the  simplified  moving  finite  difference  (MFD)  scheme 
of  Kansa  at  al.  [2].  In  that  acheme,  the  penalty  func¬ 
tion  approach  for  controlling  grid  motion  was  dropped  in 
favor  of  a  time  step  control  aod  a  static  regridding 
schema.  The  basic  modification  in  this  paper  is  that  the 
grid  velocitiea  are  part  of  an  explicit  time-integration 
schema. 

The  conservation  equations  in  a  moving  frame  have  the 
following  form  (relative  to  the  fixed  frame): 


Ht  ■"  (f*  -  vl  ♦  (Gy  -  v2  -  o, 
dx/dt  «  vl,  dy/dt  ■  v2, 


(4) 


where  the  grid  velocity,  v,  ia  arbitrary.  In  the  absence 
of  source  terms,  each  PDE  has  a  grid  velocity  for  which 
a  specific  component  of  li  is  stationary  which  is  chosen 
by 

vl  -  F^/U^  if  IUxl  *  0;  v2  «  Cy/Uy  if  I  Uyl  ?  0.  (J) 

In  special  regions  such  as  shocks,  contact  discontinui¬ 
ties^  and  rarefaction  waves,  all  the  component  grid 
velocities  are  equal.  For  the  convenience  of  dealing 
with  a  common  grid,  a  common  grid  velocity  was  deter¬ 
mined  for  the  general  ease  by  requiring  the  entire  set 
of  dependant  variables  at  a  given  location  be  atat ionary 
ia  the  least-squares  sense,  (2] .  The  time-integration 
schema,  Eq.  (3),  is  the  same,  except  that  the  A-matrix 
is  now  replaced  by  A  -  vl,  and  F*  is  replaced  by  F,  -  vU*. 
In  two  dimensions,  a  rotation  ia  performed  which  maximises 
gradients  along  the  prime  coordinate,  x'.  Then  a  trans¬ 
lation  frame  is  sought  in  which  U  is  stationary  using  Eq. 
(5). 

Continuoua  Approximations  of  Discontinuities 


Korn  and  Korn  [JJ  have  presented  several  approximate 
functions  for  Che  step  function.  The  choice  used  in  this 
paper  approximates  a  step-function  by  the  arcCan  function, 

B(x)  " 'j  (tf*  *g)  ♦  (tr“  fg)/*arctan(s(x-xo))/n 

and  s  ■  k/(xr  -  xg)  for  x  in  [xg,  x-1  and  where  fr  and 
Bg  «ta  the  endpoint  values  of  the  discontinuity,  x„  is 
the  midpoint  of  the  discontinuity,  and  s  is  the  spreading 
paranetar.  Hera  the  subscripts  r  and  t  refer  to  the  right 
(ahead)  and  laft  (behind)  of  a  discontinuity  moving  to  the 
right.  In  the  limit  of  s  going  to  infinity  as  (xr  -  xg) 
goes  to  aero,  f(x)  becomes  a  true  step  function. 

Three  mash  points  are  sufficient  to  represent  a 
"discontinuity”:  the  two  endpoints  and  the  midpoint. 
Successive  differentation  of  Eq.  (6. a)  yields  the 
following  midpoint  derivatives: 

fix,,)  •  (fr-fg)  s/X  ,f"(x„)  ■  0  .  (6b) 

At  the  endpoints,  the  rules  of  calculus  for  piecewise 
continuous  functions  (see  Korn  and  Korn  (31)  are  used  to 
determine  the  derivatives  at  cither  side  of  a 
discontinuity. 
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f'Uj)  -  f'(xj); 

f"(x4)  -  f'(x^) 

(6c) 

f(x  )  •  f*(x*); 
r  r 

f"(x  )  •  f"(x*)  * 

r  r 

(6d) 

Mot*  that  at  Cha  «nd point*  of  a  rara fact  ion  fan,  tha 
function  i*  continuoua,  but  tha  darivativas  ara  not.  Eq. 
(6c)  and  (6d)  ara  u*ad  at  tha  appropriate  location*. 

Each  component  of  F  and  U  ia  fitted  by  ita  own  atap 
function  at  a  "discontinuity™.  The  grid  velocity  ia 
found  by  uaing  the  appropriate  flux  and  dependant  vari¬ 
able  derivative*,  Eq.  (6b)  in  Eq.  (5).  The  entire  "die- 
continuity"  from  xi  to  xr  novea  at  the  grid  velo¬ 
city.  Since  one  can  atore  Che  appropriate  derivative*, 
the  aaaM  time-marching  acheae,  Eq.  (3),  can  be  uaed 
everywhere  except  at  the  boundary.  Nuawrical  experimen¬ 
tation  ha*  ahown  that  in  a  Moving  fraae,  the  firat  and 
second  order  tern*  of  Eq.  (3)  are  identically  aero  where 
the  step  function  approximation  is  used.  Using  Eq.  (6) 
in  Eq.  (3)  and  Eq.  (4)  requires  no  artificial  dissipa¬ 
tion.  Shock  ficting  is  automatically  and  correctly 
handled  by  chis  choice  of  moving  grid  velocities. 

Because  "discontinuities'*,  corner*  of  rarefaction 
fans,  and  boundaries  are  flagged,  any  shock  flow  can  be 
partitioned  into  several  distinct  piecewise  continuous 
subregions  bounded  by  flagged  and  or  boundary  points. 

Grid  motion  in  the  subregions  can  be  controlled  by 
equipartitioning  the  component  averaged  third  derivative, 
maintaining  minimum  spacing  and  well  ordering  of  paints, 
see  ( 2 J .  The  remapping  is  totally  conservative  becauae 
the  integral  of  the  old  partition  is  forced  to  equal  the 
integral  of  the  new  partition.  In  tome  circumstances, 
new  discontinuities  tend  to  develop. 

The  nonlinear  convective  terms  of  the  governing  POE't 
permit  steepening  and  eventual  discontinuity  formation. 

In  a  moving  grid  schema,  chis  phenomenon  has  bean  docu¬ 
mented  (tee  (2  and  4)).  Given  the  current  time  step, 
grid  positions,  and  velocities,  a  search  is  made  within 
a  piecewise  continuous  region  to  determine  whether  grid 
point*  would  cone  within  a  separation  of  10**  or  closer 
in  the  next  time  cycle,  and  if  gradients  are  becoming 
very  large.  If  so,  a  "discontinuity"  defined  by  three 
points,  see  Eq.  (6),  is  inserted  into  Che  computational 
domain,  and  another  flag  is  set. 

In  chis  paper,  a  nonlinear  1-D  Kiemann  solver  is  used 
co  specify  che  physics  of  wave  collisions.  The  Eiamann 
solver  is  not  used  in  any  variations  of  the  Gudunov 
scheme  co  time  advance  the  solution.  Therefore,  it  is 
instruccive  co  outline  che  use  of  the  Eiamann  solver  in 
specifying  the  resulting  physics  after  waves  collide, 
such  as  shocks  with  rigid  walls,  geometric  symaetry 
points,  shocks  with  other  shocks  or  contact  surfaces, 
etc. 

The  Riemann  solver  will  produce  one  of  five  possible 
middle  state  configurations  depending  upon  input  left 
and  right  state  densities,  velocities  and  pressures, 
located  at  Xg  and  Xr.  They  are:  (1)  middle  vacuum 
state,  (2)  a  left  and  right  rarefaccion,  (3)  a  left  shock 
and  right  rarefaction,  (4)  a  left  and  right  shock,  and 
( S)  a  left  rarefaction  and  a  right  shock.  If  the  middle 
left  and  middle  right  densities  are  unequal  then  a  den¬ 
sity  discontinuicy  exists.  Discontinuities  and  rarefac¬ 
tions  are  flagged. 

Before  advancing  the  solution  and  grid  position  co  a 
new  time,  a  search  is  made  of  all  neighboring  "discon¬ 
tinuities"  using  che  relative  separation*  and  velocities 
to  determine  whether  a  wav*  collision  is  possible  during 
the  current  time  interval.  If  some  collisions  are 
possible,  Chen  the  minimum  collision  time  interval  is 
found,  and  as  a  precaution,  it  is  multiplied  by  0.3. 

If  che  separation  of  two  existing  waves  is  less  Chan 
or  equal  to  e,  where  e  is  arbitrarily  chosen  to  be  10*’ , 
then  the  time  interval  to  complete  collision,  &ti,  can  be 
calculated.  Then  the  Riemann  solver  yields  new  states 
with  new  wav*  speeds.  New  discontinuities  and  rarefac¬ 
tions  are  flagged.  The  cim*  interval,  &t2,  to  the  same 
separation  can  be  calculated  as  well  as  the  location 
within  that  separation  of  all  the  intermediate  waves 
using  the  new  wave  velocities.  Rather  than  fixing 


the  total  number  of  grid  points,  the  number  of  grid 
points  can  be  dynamically  increased  or  decreased  after 
eh*  post  collision  state*  are  determined. 

Presentation  of  Result* 


Six  example  problems  are  presented.  These  are:  the 
Riemann  shock  Cub*  problem  in  cartesian  geometry,  the  von 
Neumann  spherical  blast  wave  problem  (31,  the  spherical 
blase  wave  problem  of  Noh  [6],  and  the  shock  tub*  pro¬ 
blem  involving  an  interacting  shock  with  a  step-down  and 
a  step-up  contact  discontinuity.  The  last  example  is  an 
infinite  planar  cwo-dimansional  shock  propagating 
obliquely  to  the  coordinate  axes. 

411  example  calculation*  were  solved  in  dimensionles* 
form.  The  solutions  have  multiple  time  and  length  scales 
which  evolve  in  time.  The  reference  length  scale  is  the 
distance  from  the  origin  at  the  left  to  the  far  right 
boundary,  scaled  to  unity.  The  reference  time  scale  is 
the  reference  length  scale  divided  by  the  initial  maximum 
speed.  Relative  to  the  initial  length  scale,  the  shock 
half  width  is  arbitrarily  chosen  to  be  10** . 

The  first  case  to  be  considered  is  the  Riemann  shock 
tube  problem  which  was  previously  calculated  by  Sod  [7] 
and  Harten  and  Hyman  [8]  with  the  initial  discontinuity 
at  x  •  1/2.  At  t  ■  0,  a  Riemann  solver  was  used  to 
initialise  the  solution  at  a  time  step  At,  on  the  do¬ 
main  (0,1).  A  rarefaction  wav*  propagates  backwards 
while  a  shock,  followed  by  a  contact  surface,  propagates 
forward.  This  system  has  an  analytic  self-similar  solu¬ 
tion.  The  shock  and  contact  surface  propagate  forwerd  at 
conatant  velocities.  The  rarefaction  fan  propagates 
backwards  along  the  local  u-a  characteristic,  where  u  and 
a  are  eh*  local  particle  and  sound  velocities. 

Figure  1  show*  the  numerical  solution  for  the  density 
of  y  *  1.4  gas  at  e  ■  0.23.  The  diamonds  in  these 
figure*  and  all  subsequent  figures  represent  the  actual 
solution  at  tha  indicated  nodal  position.  Not*  that  be¬ 
cause  the  conservative  variables  (see  Eq.  2)  were  used, 
the  gas  velocity  and  pressure  at  a  discontinuity  may  not 
necessarily  be  the  midpoint  values. 

Tha  second  problem  attempted  was  the  von  Neumann  [3] 
infinite  shock,  spherical  blast  wave  problem.  Ahead  of 
the  shock,  p  •  u  “  0  and  0*1.  The  pressure  behind 
tha  shock  at  radius,  r  *  l,  was  specified  to  be  100.  A 
special  routine  was  used  to  initialize  the  solution  be¬ 
hind  the  shock. 

From  the  analysis  of  Zal'dovich  et  al.  [9)  and 
Thompson  [10],  the  following  relations  can  be  obtained 
for  ay*  5/3  gas  in  terms  of  new  veriables: 

0  *  4  p,  u  *  (r/t)  0,  p  *  (r/t)2p.  (7) 

Using  the  above  relations  in  a  and  E  as  well  as  the 
fluxes,  the  following  relation  was  observed: 

Ut  ♦  (l/t)(l/r2)(r2  F)r  -  0.  (8) 

Because  of  Che  (1/t)  multiplying  che  flux  transport  terms, 
the  integrating  factor  of  Eq.  (3)  is  En(t). 

Tha  results  for  pressure  are  shown  in  Fig.  2  at  cwo 
different  times.  In  the  moving  frame,  only  density  is 
self-similar.  The  velocity  and  pressure  profiles  decay. 
The  in(c)  integrating  factor  gives  the  correct  temporal 
behavior  without  unnecesaarily  small  time  step*. 

The  next  problem  co  be  discussed  is  Noh's  [6]  spheri¬ 
cally  divergent  shock  problem.  At  time  t  •  0,  a  y  *  5/3 
gas  sc  zero  pressure,  density  of  one,  converges  within 
a  sphere  at  a  gas  velocity  of  -1.  Then  at  the  center,  a 
shock  develops  with  a  density  of  64,  gas  velocity  of  zero, 
and  pressure  of  64/3  moving  outwards  at  a  shock  velocity 
of  1/3.  Noh's  challenge  is  to  run  the  problem  co  a  time 
of  t  *  0.6. 

In  order  to  solve  the  shock  scace  due  co  the  focusing 
of  gas  at  the  origin,  the  Riemann  solver  must  be  modified 
for  infinite  shocks  and  cylindrical  and  spherical  geome¬ 
try.  As  soon  as  Che  gat  moves  inward,  it  is  adiabatically 
compressed  by  the  spherical  focusing.  Since  the  inicial 
pressure  it  zero,  the  subsequent  pressure  profile  inward, 
but  not  at  the  origin  it  alto  zero.  The  inward  gas  motion 


it  self-similar  with  a  constant  gat  velocity  of  u  *  -1. 
and  the  matt  conservation  aquation  can  ba  integrated  Co 
yield 

t>  «  (1  ♦  t/r)J  (8) 

for  r>o 

with  eha  boundary  condicion 

r  ■  -t,  at  J  ■  1  . 

Becaute  of  tpherical  symmecry,  the  gat  valocicy  mute 
ba  aero  at  the  origin.  Uting  the  notation  that  the  sub- 
tcriptt  t  and  r  refer  to  the  left  and  right  ttatea  ahead 
of  the  shock,  and  •  to  Che  new  middle  shock  ttace,  the 
known  quantitiet  are  un,  Pt  "  Pr  "  ur  ”  “1  ■  ug  ”1  • 

The  unknown  quantitiet  are  Pa>  pa,  uSr  •  -uSg,  and  pr  ■  Pg. 
Following  Couranc  et  al.  [11J,  one  finde  for  a  y  1  5/3  gat 
that  pr  •  16,  Pa  *  6*,  utf  ”  1/3,  and  pa  •  64/3.  Fig.  3 
is  the  density  calculated  to  t  •  0.6.  Note  that  the 
densicy  has  developed  a  new  cootact  surface  at  the  tail 
before  it  abruptly  drops  to  unity.  This  discontinuity 
was  inserted  as  discussed  previously. 

The  next  problea  in  slab  geoaetry  is  adapted  froa 
Harlow  and  Aasden  [12).  At  c  *  0,  an  incooing  gat  of 
velocity  u  “  -1,  p  •  0  hits  a  rigid  wall  located  at  x  *  0- 
The  inicial  gat  density  is  2  for  x  1/2  and  unity  for 
1/2  <  x  <_  l.  The  Riemann  solver  produces  an  outgoing 
shock  noving  at  a  velocity  of  1/3  and  an  incoaing  contact 
discontinuity  aoving  at  a  velocity  of  -1.  Figure  4  shows 
the  density,  at  c  ■  0.2  after  the  shock  formation.  Figure 
4  alto  shows  the  density  at  e  •  0.3169  jutt  prior  to  the 
shock  colliding  with  the  contact  discontinuity. 

At  the  time  of  collision,  the  Kicaann  solver  is  used 
again,  yielding  a  left  traveling  rarefaction,  and  a  right 
traveling  contact  discontinuity  and  a  shock  with  veloci¬ 
ties  0.166  and  0.3347,  respectively.  The  post  collision 
states  were  resolved  adequately  by  adaptive  mesh  refine¬ 
ment.  Figure  4  also  shows  the  density  at  t  *  0.4184. 

The  fifth  problea  to  be  considered  is  also  taken  froa 
Harlow  and  Aasden  [12).  At  t  •  0,  the  initial  state  was 
defined  to  have  a  pressure  equal  to  lO'*’ ,  an  incoaing 
gas  velocity  of  -1,  and  a  denaity  of  1  for  x  <_  1/2  and 
8/3  for  x  <  1/2  <  1.  At  x  ■  0,  there  is  a  rigid  wall. 

A  {low  with  a  denaity  jwp  hits  a  rigid  wall,  and  a  shock 
is  formed.  The  shock  proceeds  to  the  right  and  collides 
with  a  step-up  concact  discontinuity.  The  Rieaann  solver, 
at  the  time  of  collision,  gives  as  solutions  a  right  and 
left  traveling  shock.  Figure  3  shows  tha  shock  traveling 
at  t  •  0.2  to  the  right  with  a  valocity  Vg  -  0.333  and 
a  contact  surface  moving  to  Che  left  at  the  valocity  of 
-1.  Figure  3  also  shows  the  new  states  after  collision 
at  t  ■  0.31867.  Although  this  appears  to  be  an  overshoot, 
the  new  states  are  spread  over  a  distance  of  0.001.  The 
new  state  has  two  shocks,  the  left  shock  has  a  velocity 
of  -0.909  and  the  right  shock  has  a  velocity  of  0.0336. 

The  middle  contact  has  a  velocity  of  -0.223. 

Figure  6  shows  Che  densicy  at  t  •  0.4186  with  left 
shock  approaching  Che  wall.  Ac  t  ■  0.463,  Fig.  6  shows 
the  densicy  plot  after  wall  collision.  The  shock  hits 
Che  wall  proceeding  forward  with  a  velocity  of  0.7602. 
Figure  7  at  t  ■  0.363  shows  the  density  profile  just 
prior  to  the  collision  of  the  contact  with  the  lefc  shock 
which  reflected  off  Che  well.  Figure  7  also  shows  the 
profile  after  Che  shock  from  the  wall  collided  with  the 
concact  at  t  *  0.703. 

These  calculations  were  extended  at  will  without  any 
evidence  of  difficulty.  The  celculations  were  arbitrarily 
stopped  at  t  •  0.803.  No  instabilities  heve  occurred 
because  accurate  solutions  from  the  tiemann  solver  were 
used  to  define  the  collision  states,  a  moving  grid  scheme 
was  used,  and  a  nathematically-consistent  approximate 
step  function  was  used  without  artificial  viscosity.  In 
this  particular  calculation,  no  CFL  time  restrictions 
were  necessary  since  Che  first  end  second  order  tine 
corrections  (see  Eq.  9.b)  were  identically  sero.  The 
time  scale  estimate  used  is  t  ■  HUH/llV*F  -  v*VUU. 

When  Che  denominator  vanishes,  Che  time  scale  is  arbi¬ 
trarily  large,  however,  the  usual  CFL  condition  is 
recovered  when  the  grid  mocion  is  turned  off. 


The  last  problem  is  a  demonstration  thaC  this  noving 
grid  scheme  does  indeed  extend  to  higher  dimensions.  On 
a  unit  square,  an  infinite  shock  propagates  obliquely  to  I 
tha  coordinate  axes.  To  the  left  of  the  shock  for  a 
y  -  5/3  ideal  gas,  p  -  4/3,  p  »  4,  u  •  v  -  0;  to  the 
right  of  the  shock,  p  ■  0,  P  *  1,  u  •  -l,  v  •  -  5/8. 

The  initial  pressure  profile  is  shown  in  Figure  8. 

Two  grids  were  used.  An  underlying  coarse  grid  was 
used  with  Ax  ■  1/6  and  Ay  "  1/5*  The  shock  region 
is  defined  on  a  mesh  of  length  7-11  grid  points  in  the 
tangential  direction  and  a  width  of  3  mesh  points  of 
total  extent  10“4  in  the  notmal  direction.  A  two 
dimensional  extension  of  a  second  order  time  marching 
scheme  was  used  throughout  the  domain,  except  at  the 
boundaries  which  were  open.  At  the  shock,  a  rotation  was 
performed  into  the  normal  direction,  and  then  a  tranala- 
tional  frame  was  found  in  which  the  conservative  depen¬ 
dent  variables  were  steady.  The  pressure  profiles  at  t 
>  0  and  t  ■  0.56  are  shown.  Note  that  the  shock  remains 
sharp.  These  calculations  were  run  at  a  CFL  number  of 
200.  Figure  9  shows  these  results. 

The  future  plans  are  to  incorporate  the  two  dimen¬ 
sional  shock  interaction  analysis  of  Courant  [11]. 

Because  of  the  network  of  shock  interactions,  i.e., 
triple  points,  bubbles,  etc.  expected  to  arise  from  two 
dimensional  interacting  shocks,  a  rectangular  mesh  may 
become  too  cumbersome  and  scattered  data  interpolation 
schemes  might  be  most  useful. 
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